/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#ifndef UHM_INTERF_PARDISO_HXX
#define UHM_INTERF_PARDISO_HXX

namespace uhm {

  namespace interf {
    typedef class Sparse_*  Sparse;
    typedef class Pardiso_* Pardiso;

    // ----------------------------------------------------------------
    // ** PARDISO
    class Pardiso_ {
    private:
    protected:

      int datatype, show_n_rhs;

      // pardiso control parameter

      // execution step
      // 11 - analysis
      // 12 - analysis + factorization
      // 13 - analysis + factorization + solve + iterative refinement
      // 22 - factorization
      // 23 - factorization + solve + iterative refinement
      // 33 - solve + iterative refinement
      // 0  - release internal memory for LU matrix number MNUM
      // -1 - release all internal memory for all matrices
      int phase;

      // internal data address pointer never touch
      void  *pt[64];

      //  1 - real structrue sym
      //  2 - real sym pos def
      // -2 - real sym indef
      //  3 - complex structrue sym
      //  4 - complex hermitian pos def
      // -4 - complex and hermitian indef
      //  6 - complex and sym
      // 11 - real and nonsym
      // 13 - complex and nonsym
      int    mtype;      

      //  0 - direct , 1 - multi recursive iterative solver
      int    solver;

      // parameters - see manual
      int    iparm[64];
      double dparm[64];
    
      int    maxfct, mnum, msglvl;

      int    num_procs;

      // matrix CSR format
      int n_dof, n_rhs, n_nz;
      std::vector<int> ia, ja;
      std::vector<double> a, b, x;

      void _init(int datatype) {
        assert(datatype == UHM_REAL ||
               datatype == UHM_COMPLEX);
        this->datatype = datatype;
        this->show_n_rhs = 0;
      }
    
    public:

      // methods
      Pardiso_() { _init(UHM_REAL); }
      Pardiso_(int datatype) { _init(datatype); }
      virtual ~Pardiso_() { }

      bool export_matrix(Mesh m);

      void disp();
      void disp( FILE *stream );


      bool is_complex();
      void set_show_n_rhs(int show_n_rhs);
      void set_env(int mtype, int solver, int maxfct, int mnum, int msglvl);
      void set_phase(int phase);
      void set_iparm(int idx, int val);
      void set_dparm(int idx, double val);
      void set_sparse_matrix(Sparse sp, int is_sym);

      int    get_iparm(int idx);
      double get_dparm(int idx);

      int  get_n_dof();    
      int  get_n_rhs();    
      int  get_n_nonzero();

      bool validate_sparse_matrix();
      int  run();

      bool init();
      bool analyze();
      bool decompose();
      bool solve();
      bool finalize();
    };

    // ----------------------------------------------------------------
  }
}
#endif
